function M = mprecond(f,x)
% function M = mprecond(f,x)
%   this function produces a matrix M for use as a preconditioner in 
%   solving Newton step, i.e.
%          dx = -inv(J)*M*F(x)
%   It is based on the paper by J.L. Hueso et al.
n   = length(x);
M   = zeros(n);
del = 1e-5;

Fo  = feval(f,x);
Jo  = diffjac(x,f,Fo);%jacob(f,x,Fo);
JFo = Jo\Fo;

for i = 1:n
    w       = zeros(n,1);
    w(i)    = del;
    F       = feval(f,x+w);
    J       = jacob(f,x+w,F);
    JF      = J\F;
    M(i,i)  = del/(JF(i)-JFo(i));
    if (M(i,i) < 1.0)
        M(i,i) = 1.0;
    end
end


